Introduction
This file contains essential commands from Chapter 7: Exploratory data analysis of the textbook r4ds and corresponding examples and exercises. A command is considered “essential” when you really need to know it and need to know how to use it to succeed in this course.
All ds4psy essentials so far:
| Nr. | Topic |
|---|---|
| 0. | Syllabus |
| 1. | Basic R concepts and commands |
| 2. | Visualizing data |
| 3. | Transforming data |
| 4. | Exploring data (EDA) |
| +. | Datasets |
Course coordinates
- Taught at the University of Konstanz by Hansjörg Neth (h.neth@uni.kn, SPDS, office D507).
- Winter 2018/2019: Mondays, 13:30–15:00, C511.
- Links to current course syllabus | ZeUS | Ilias
Preparations
Create an R script (.R) or an R-Markdown file (.Rmd) and load the R packages of the tidyverse. (Hint: Structure your script by inserting spaces, meaningful comments, and sections.)
## Exploring data (EDA) | ds4psy
## 2018 11 26
## ----------------------------
## Preparations: ----------
library(tidyverse)
## 1. Topic: ----------
# etc.
## End of file (eof). ---------- Exploring data
Introduction
In the following, we refine and use what we have learned so far to explore data. Practically, this session combines what we have learned about ggplot2 (in Chapter 3: Data visualization) and about dplyr (in Chapter 5: Data transformation) to explore datasets.
Key concepts
Important concepts in this session include:
- missing values (
NA) vs. unusual and extreme values (e.g., outliers), - different types of variables (e.g., continuous vs. categorical variables),
- relationships between variables (of different types),
- trends (i.e., developments over time), and
- many different types of plots.
See Chapter 7: Exploratory data analysis (EDA) and the links provided below for additional information.
What is EDA?
According to Grolemund & Wickham (2017, Chapter 7) exploratory data analysis (EDA) primarily is a state of mind. Rather than following a strict set of rules, EDA is an initial phase of familiarising yourself with a new data set by actively engaging in an iterative cycle during which you
- Generate questions about data.
- Search for answers by visualising, transforming, and modelling data.
- Use answers to refine your questions and/or generate new questions.
Philosophically speaking, EDA is a data scientist’s way of doing hermeneutics (see Wikipedia or The Stanford Encyclopedia of Philosophy for definitions) to get a grip on data.
The goal and purpose of EDA is to gain an overview of a dataset. This includes an idea of its dimensions and the number and types of variables contained in the data, but also a more detailed idea of the distribution of variables, their potential relations to each other, and potential problems (e.g., missing values or outliers).
When using the tools provided by the tidyverse, the fastest way to gain insights into a dataset is a series of dplyr calls and ggplot2 graphs. However, creating good graphs is both an art and a craft. The key to creating good graphs requires answering 2 sets of questions:
Knowing the intended type of plot. This includes answering functional questions like
- What is the goal or purpose of this plot?
- What are possible plot types for this purpose?
- Which of these would be the most appropriate plot here?
Knowing the number and type of variables to be plotted. This includes answering data-related questions like
- How many variables are to be plotted and how are they mapped to dimensions and aesthetics?
- Are these variables categorical or continuous?
- Do some variables control or qualify (e.g., group) the values of others?
- How many variables are to be plotted and how are they mapped to dimensions and aesthetics?
Even when these questions are answered, creating beautiful and informative graphs with ggplot requires dedicated practice, experience, and trial-and-error experimentation. In addition, an new dataset is rarely in a condition and shape to directly allow plotting. Instead, we typically have to interleave commands for plotting and transforming data to wrangle variables into specific shapes or types. Thus, calls to ggplot2 usually occur in combination with other tidyverse commands (stemming from dplyr, forcats, tidyr, readr, tibble, etc.).
What needs to be done in any specific case depends on the details of the data and your current goals. Also, there is no pre-defined end to an EDA and no clear boundary between the processes of exploration and the confirmation (or falsification) of expectations. Typically, social scientists use EDA to check and understand a dataset, before using statistics to test specific hypotheses.
While any actual EDA is tailored to specific features of the data and current research goals, the following sections highlight some common themes that occur in most cases. As we will see, the packages dplyr and ggplot provide a powerful set of tools that allow us to ask and answer questions about data. Ideally, you will soon experience that actively engaging in EDA really feels like “doing research”, as it requires us to probe and scrutinize data and follow up on the ideas and hypotheses that we may discover along the way.
To illustrate the process and spirit of EDA, we will use a real and fairly complex data set, which was collected to measure the short- and long-term effects of positive psychology interventions (see Dataset: Positive Psychology for details).1 And while it is understandable that you’re mostly focusing on the code and commands, let’s try to stay curious about the results that we may find.
Principles and practices
In the following, we will explain some principles that endorse and promote the ideal of transparent data analysis and reproducible research. While such practices are indispensable when working in a team of colleagues and the wider scientific community, organizing your workflow in a more consistent fashion is also beneficial for your other projects and your future self.
Setting the stage
Before embarking on any actual analysis, we should pay tribute to 2 principles that seem so simple that it’s easy to overlook their benefits:
- Principle 1: Structure and comment your analysis.
While adhering to this principle may seem superfluous or trivial at first, it is crucial as analyses get longer and more complicated (e.g., distributed over multiple datasets and scripts). A consistent document structure, transparent object names, and clear and informative comments are indispensable when sharing your data or scripts with others. Good comments often structure a longer document into parts and briefly explain the content of each part (what), but will primarily focus on the goals and reasons for your decisions (i.e., explain why, rather than how you did something).
- Principle 2: Start with a clean slate and explicitly load all data and all packages required in this analysis.
Although RStudio IDE provides options for saving your workspace and for loading data or packages by clicking buttons or checking boxes, doing so renders the process of your analysis intransparent for anyone not observing all your actions. To make sure that others and yourself can repeat the same sequence of steps tomorrow or next year, it is advisable to always start with a clean slate and explicitly state which data and packages are required for the current analysis (unless there are good reasons to deviate from them):2
Cleaning up
Clean your workspace and define some custom objects (like colors or functions):
# Housekeeping:
rm(list = ls()) # cleans ALL objects in current R environment (without asking for confirmation)!
# Customizations:
seeblau <- rgb(0, 169, 224, names = "seeblau", maxColorValue = 255) # seeblau.4 of uni.kn color scheme
# source(file = "my_custom_functions.R")Loading packages and data
Load all required packages and data files (see Dataset: Positive Psychology for details):
# Load packages:
library(tidyverse)
library(knitr)
# Load csv-data files (from online links):
# See <http://rpository.com/ds4psy/essentials/datasets.html> for details.
# 1. Participant data:
posPsy_p_info <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_participants.csv")
# 2. Original DVs in long format:
AHI_CESD <- read_csv(file = "http://rpository.com/ds4psy/data/posPsy_AHI_CESD.csv")
# 4. Corrected version of all data in wide format:
posPsy_wide <- readr::read_csv(file = "http://rpository.com/ds4psy/data/posPsy_data_wide.csv")- Principle 3: Make copies (and copies of copies) of your data.
It is only human to make occasonal errors. To make sure that errors are not too costly, a good habit is to occasionally save intermediate steps. This refers to both your data files and your scripts for analyzing them. Although it may be a good idea to always keep a ``current master version’’ of a data file, it is also advisable to occasionally copy your current data and then work on the copy (e.g., before trying out some unfamiliar analysis or command). In R, creating copies and then working on them is very easy — and can even save yourself from repeatedly typing complicated object names. Most importantly, when working on a copy of your data, you can always recover the last sound version if things go terribly wrong:
df <- posPsy_wide # copy data (and simplify name)
dim(df) # 295 cases x 294 variables
#> [1] 295 294
sum(is.na(df)) # 37440 missing values
#> [1] 37440
df$new <- NA # create a new column with NA values
df$new # check the new column
#> [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [24] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [47] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [70] NA NA NA NA NA NA
#> [ reached getOption("max.print") -- omitted 220 entries ]
df <- NA # create a new column with NA values
# Ooops...
df # looks like we accidentally killed our data!
#> [1] NA
# But do not despair:
df <- posPsy_wide # Here it is again:
dim(df) # 295 cases x 294 variables
#> [1] 295 294During a long and complicated analysis, it is advisable to save an external copy of your current data when having reached an important intermediate step. While there are many different formats in which data files can be stored, a good option for most files is csv (comma-separated-values), which can easily be read by most humans and machines. Importantly, always verify that the external copy preserves the key information contained in your current original:
# Write out current data (in csv-format):
write_csv(df, path = "my_precious_data.csv")
# Re-read external data (into a different object):
df_2 <- read_csv(file = "my_precious_data.csv")
# Verify that original and copy are identical:
all.equal(df, df_2)
#> [1] TRUEScreening data
Screening data involves checking the dimensions of data, the types of variables, missing or atypical values, and the validity of observations.
Basic properties
What do you want to know immediately after loading a data file?
- Principle 4: Know your data (variables and observations).
This principle is a short version of the following: You really want to know
- the dimensions of your data,
- the types of your variables (columns), and
- the semantics of your observations (rows). Here, we first inspect the original data file, which we read in from http://rpository.com/ds4psy/data/posPsy_AHI_CESD.csv above:
df <- AHI_CESD # copy the data
dim(df) # 992 cases x 50 variables
#> [1] 992 50
df <- as_tibble(df) # (in case df isn't a tibble already)
# Get an initial overview:
df
#> # A tibble: 992 x 50
#> id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05
#> <int> <int> <dbl> <int> <int> <int> <int> <int> <int>
#> 1 1 0 0 4 2 3 2 3 3
#> 2 1 1 11.8 4 3 3 4 3 3
#> 3 2 0 0 1 3 4 3 4 2
#> 4 2 1 8.02 1 3 4 4 4 3
#> 5 2 2 14.3 1 3 4 4 4 3
#> 6 2 3 32.0 1 3 4 4 4 4
#> 7 2 4 92.2 1 3 3 2 3 3
#> 8 2 5 182. 1 3 3 3 4 2
#> 9 3 0 0 4 3 3 2 4 2
#> 10 3 2 16.4 4 3 3 3 4 4
#> # ... with 982 more rows, and 41 more variables: ahi06 <int>, ahi07 <int>,
#> # ahi08 <int>, ahi09 <int>, ahi10 <int>, ahi11 <int>, ahi12 <int>,
#> # ahi13 <int>, ahi14 <int>, ahi15 <int>, ahi16 <int>, ahi17 <int>,
#> # ahi18 <int>, ahi19 <int>, ahi20 <int>, ahi21 <int>, ahi22 <int>,
#> # ahi23 <int>, ahi24 <int>, cesd01 <int>, cesd02 <int>, cesd03 <int>,
#> # cesd04 <int>, cesd05 <int>, cesd06 <int>, cesd07 <int>, cesd08 <int>,
#> # cesd09 <int>, cesd10 <int>, cesd11 <int>, cesd12 <int>, cesd13 <int>,
#> # cesd14 <int>, cesd15 <int>, cesd16 <int>, cesd17 <int>, cesd18 <int>,
#> # cesd19 <int>, cesd20 <int>, ahiTotal <int>, cesdTotal <int>
## Other ways to probe df:
# names(df) # Note variable names
# glimpse(df) # Note types of variables
# summary(df) # Note range of valuesWe note that AHI_CESD contains 992 observations (rows) and 50 variables (columns). Most of the variables are integers and – judging from their names – many belong together (in blocks). Their names and ranges suggest that they stem from questionnaires or scales with fixed answer categories (e.g., from 1 to 5 for ahi__, and from 1 to 4 for cesd__). Overall, a row of data is a measurement of one participant (characterised by its id and an intervention value) at one occasion (characterised by the value of occasion and elapsed.days) on two scales (ahi and cesd).
Unusual values
- Principle 5: Know and deal with unusual values.
Unusual values include:
- missing values (`NA` or `-77`, etc.),
- extreme values (outliers),
- other unusual values (e.g., unexpected and impossible values). Missing values
What are missing values? In R, missing values are identified by NA. Note that NA is different from NULL: NULL represents the null object (i.e., something is undefined). By contrast, NA means that some value is absent.
Both NA and NULL are yet to be distinguished from NaN values.
## Checking for NA, NULL, and NaN:
# NA:
is.na(NA) # TRUE
#> [1] TRUE
is.na(NULL) # 0 (NULL)
#> logical(0)
is.na(NaN) # TRUE!
#> [1] TRUE
# NULL:
is.null(NULL) # TRUE
#> [1] TRUE
is.null(NA) # FALSE
#> [1] FALSE
is.null(NaN) # FALSE
#> [1] FALSE
# NaN:
is.nan(NaN) # TRUE
#> [1] TRUE
is.nan(0/0) # TRUE
#> [1] TRUE
is.nan(1/0) # FALSE, as it is +Inf
#> [1] FALSE
is.nan(-1/0) # FALSE, as it is -Inf
#> [1] FALSE
is.nan(0/1) # FALSE (as it is 0)
#> [1] FALSE
is.nan(NA) # FALSE
#> [1] FALSE
is.nan(NULL) # 0 (logical)
#> logical(0)
# Note different modes:
all.equal(NULL, NA) # Note: NULL vs. logical
#> [1] "target is NULL, current is logical"
all.equal("", NA) # Note: character vs. logical
#> [1] "Modes: character, logical"
#> [2] "target is character, current is logical"
all.equal(NA, NaN) # Note: logical vs. numeric
#> [1] "Modes: logical, numeric"
#> [2] "target is logical, current is numeric"Missing data values require attention and special treatment. They should be identified, counted and/or visualized, and often removed or recoded (e.g., replaced by other values).
Counting NA values:
df <- AHI_CESD # copy the data
is.na(df) # asks is.na() for every value in df
#> id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> ahi05 ahi06 ahi07 ahi08 ahi09 ahi10 ahi11 ahi12 ahi13 ahi14 ahi15
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> ahi16 ahi17 ahi18 ahi19 ahi20 ahi21 ahi22 ahi23 ahi24 cesd01 cesd02
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> cesd03 cesd04 cesd05 cesd06 cesd07 cesd08 cesd09 cesd10 cesd11
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> cesd12 cesd13 cesd14 cesd15 cesd16 cesd17 cesd18 cesd19 cesd20
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> ahiTotal cesdTotal
#> [1,] FALSE FALSE
#> [ reached getOption("max.print") -- omitted 991 rows ]
sum(is.na(df)) # sum of all instances of TRUE
#> [1] 0Since df does not seem to include any NA values, we create a tibble tb that does:
# Create a df with NA values:
set.seed(42) # for replicability
nrows <- 6
ncols <- 6
tb <- as_tibble(matrix(sample(c(1:12, -66, -77), nrows * ncols, replace = TRUE), nrows, ncols)) # create some df
tb[tb > 9] <- NA # SET certain values to NA
# tbCounting and recoding NA values:
# Count NA values:
sum(is.na(tb)) # => 10 NA values
#> [1] 10
# The function `complete.cases(x)` returns a logical vector indicating which cases in tb are complete:
complete.cases(tb) # test every row/case for completeness
#> [1] FALSE FALSE FALSE FALSE TRUE FALSE
sum(complete.cases(tb)) # count cases/rows with complete cases
#> [1] 1
which(complete.cases(tb)) # indices of case(s)/row(s) with complete cases
#> [1] 5
tb[complete.cases(tb), ] # list all complete rows/cases in tb
#> # A tibble: 1 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 9 7 -77 -77 7 1
tb[!complete.cases(tb), ] # list all rows/cases with NA values in tb
#> # A tibble: 5 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -66 NA -77 7 2 NA
#> 2 -77 2 4 8 8 NA
#> 3 5 NA 7 -66 6 6
#> 4 NA NA -77 2 -66 NA
#> 5 8 NA 2 -77 NA NA
# Recode all instances of NA as -99:
tb
#> # A tibble: 6 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -66 NA -77 7 2 NA
#> 2 -77 2 4 8 8 NA
#> 3 5 NA 7 -66 6 6
#> 4 NA NA -77 2 -66 NA
#> 5 9 7 -77 -77 7 1
#> 6 8 NA 2 -77 NA NA
tb[is.na(tb)] <- -99 # recode NA values as - 99
tb
#> # A tibble: 6 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -66 -99 -77 7 2 -99
#> 2 -77 2 4 8 8 -99
#> 3 5 -99 7 -66 6 6
#> 4 -99 -99 -77 2 -66 -99
#> 5 9 7 -77 -77 7 1
#> 6 8 -99 2 -77 -99 -99More frequently, special values (like -66, -77 and -99) indicate missing values. To deal with them in R, we recode all instances of these values in tb as NA:
tb
#> # A tibble: 6 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -66 -99 -77 7 2 -99
#> 2 -77 2 4 8 8 -99
#> 3 5 -99 7 -66 6 6
#> 4 -99 -99 -77 2 -66 -99
#> 5 9 7 -77 -77 7 1
#> 6 8 -99 2 -77 -99 -99
# Recode -66, -77, and -99 as NA:
tb[tb == -66 | tb == -77 | tb == -99 ] <- NA
tb
#> # A tibble: 6 x 6
#> V1 V2 V3 V4 V5 V6
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 NA NA NA 7 2 NA
#> 2 NA 2 4 8 8 NA
#> 3 5 NA 7 NA 6 6
#> 4 NA NA NA 2 NA NA
#> 5 9 7 NA NA 7 1
#> 6 8 NA 2 NA NA NA
sum(is.na(tb)) # => 19 NA values
#> [1] 19For more sophisticated ways of dealing with (and visualising) NA values, see the R packages mice (Multivariate Imputation by Chained Equations), VIM (Visualization and Imputation of Missing Values) and Amelia II.
Other unusual values: Unexpected values, outliers, etc.
While dealing with missing values is a routine task, finding other unusual values involves some hard thinking, but can also be fun, as it requires the state of mind of a detective who gradually uncovers facts – and aims to ultimately reveal the truth – by questioning a witness or suspect. In principle, we can detect atypical values and outliers (see Exercise 3 of WPA03 for different definitions) by counting values, computing descriptive statistics, or by checking the distributions of raw values. To illustrate this, let’s inspect the data of AHI_CESD further. Since participants in this study were measured repeatedly over a range of several months, two good questions to start with are:
- How many participants are there for each occasion?
- How many occasions are there for each participant?
The first question concerns a trend (i.e., development over time) and can easily be answered by using dplyr or ggplot:
# How many participants are there for each occasion?
# Data:
df <- AHI_CESD
# (1) Summary table of grouped counts:
# Using dplyr pipe:
id_by_occ <- df %>%
group_by(occasion) %>%
count()
id_by_occ
#> # A tibble: 6 x 2
#> # Groups: occasion [6]
#> occasion n
#> <int> <int>
#> 1 0 295
#> 2 1 147
#> 3 2 157
#> 4 3 139
#> 5 4 134
#> 6 5 120
# (2) Plot count of IDs by occasion as bar plot:
# (a) Using raw data as input:
ggplot(df, aes(x = occasion)) +
geom_bar(fill = seeblau) +
labs(title = "Number of participants by occasion (from raw data)") +
theme_bw()
# (b) Using id_by_occ from (1) as input:
ggplot(id_by_occ, aes(x = occasion)) +
geom_bar(aes(y = n), stat = "identity", fill = seeblau) +
labs(title = "Number of participants by occasion (from summary table)") +
theme_bw()
Note the similarities and differences between these two bar plots. Although they look almost the same (except for the label on the y-axis and their title), they use completely different data as inputs. Which of the 2 plot versions would make it easier to add additional information (like error bars, a line indicating the mean counts of participants, or text labels showing the count values for each category)?
# How many occasions are there for each participant?
# Data:
df <- AHI_CESD
# Graphical solution (with ggplot):
ggplot(df) +
geom_bar(aes(x = id), fill = seeblau) +
labs(title = "Number of occasions by participant (raw)") +
theme_bw()
# This does the trick, but looks messy.
# A cleaner solution uses 2 steps:
# (1) Summary table of grouped counts:
occ_by_id <- df %>%
group_by(id) %>%
count()
occ_by_id
#> # A tibble: 295 x 2
#> # Groups: id [295]
#> id n
#> <int> <int>
#> 1 1 2
#> 2 2 6
#> 3 3 2
#> 4 4 4
#> 5 5 5
#> 6 6 3
#> 7 7 1
#> 8 8 6
#> 9 9 1
#> 10 10 1
#> # ... with 285 more rows
# (2) Plot summary table:
# (a) re-create previous plot (by now from occ_by_id):
ggplot(occ_by_id) +
geom_bar(aes(x = id, y = n), stat = "identity", fill = seeblau) +
labs(title = "Number of occasions by participant (2a)") +
theme_bw()
# (b) reordering id by count n:
ggplot(occ_by_id) +
geom_bar(aes(x = reorder(id, n), y = n), stat = "identity", fill = seeblau) +
labs(title = "Number of occasions by participant (2b)") +
theme_bw()
# (c) more efficient solution:
occ_by_id_2 <- occ_by_id %>%
arrange(n)
occ_by_id_2
#> # A tibble: 295 x 2
#> # Groups: id [295]
#> id n
#> <int> <int>
#> 1 7 1
#> 2 9 1
#> 3 10 1
#> 4 15 1
#> 5 16 1
#> 6 17 1
#> 7 20 1
#> 8 23 1
#> 9 25 1
#> 10 26 1
#> # ... with 285 more rows
## (+) Add rownames (1:n) as a column:
# occ_by_id_2 <- rownames_to_column(occ_by_id_2) # rowname is character variable!
## (+) Add a variable row that contains index of current row:
occ_by_id_2$row <- 1:nrow(occ_by_id_2)
occ_by_id_2
#> # A tibble: 295 x 3
#> # Groups: id [295]
#> id n row
#> <int> <int> <int>
#> 1 7 1 1
#> 2 9 1 2
#> 3 10 1 3
#> 4 15 1 4
#> 5 16 1 5
#> 6 17 1 6
#> 7 20 1 7
#> 8 23 1 8
#> 9 25 1 9
#> 10 26 1 10
#> # ... with 285 more rows
ggplot(occ_by_id_2) +
geom_bar(aes(x = row, y = n), stat = "identity", fill = seeblau) +
labs(title = "Number of occasions by participant (2c)") +
theme_bw()
Our exploration shows that every participant was measured on at least 1 and at most 6 occasions (check range(occ_by_id$n) to verify this). This raises two additional questions:
- Was every participant measured on occasion 0 (i.e., the pre-test)?
- Was every participant measured only once per occasion?
The first question may seem a bit picky, but do you really know that nobody showed up late (i.e., missed occasion 0) for the study? Actually, we do already know this, since we counted the number of participants and the number of participants per occasion above:
- the study sample contains 295 participants (check
nrow(posPsy_p_info)), and - the count of participants per occasion showed a value of 295 for occasion 0 in
id_by_occ.
For the record, we testify that:
nrow(posPsy_p_info) # 295 participants
#> [1] 295
id_by_occ$n[1] # 295 participants counted (as n) in 1st line of id_by_occ
#> [1] 295
nrow(posPsy_p_info) == id_by_occ$n[1] # TRUE (qed)
#> [1] TRUETo answer the second question, we can count the number of lines in df per id and occasion and then check whether any unexpected values (different from 1, i.e., n != 1) occur:
# Summary table (with dplyr):
id_occ <- df %>%
group_by(id, occasion) %>%
count()
id_occ
#> # A tibble: 990 x 3
#> # Groups: id, occasion [990]
#> id occasion n
#> <int> <int> <int>
#> 1 1 0 1
#> 2 1 1 1
#> 3 2 0 1
#> 4 2 1 1
#> 5 2 2 1
#> 6 2 3 1
#> 7 2 4 1
#> 8 2 5 1
#> 9 3 0 1
#> 10 3 2 1
#> # ... with 980 more rows
summary(id_occ$occasion)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000 0.000 2.000 2.028 4.000 5.000
summary(id_occ$n) # Max is 2 (not 1)!
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.000 1.000 1.000 1.002 1.000 2.000
# Do some occasions occur with other counts than 1?
id_occ %>%
filter(n != 1)
#> # A tibble: 2 x 3
#> # Groups: id, occasion [2]
#> id occasion n
#> <int> <int> <int>
#> 1 8 2 2
#> 2 64 4 2
# => Participants with id of 8 and 64:
# 2 instances of occasion 2 and 4, respectively.Importantly, 2 participants (8 and 64) are counted twice for an occasion (2 and 4, respectively).
Compare our counts of id_by_occ with Table 1 (p. 4, of Woodworth et al., 2018), which also shows the number of participants who responded on each of the 6 measurement occasions): As the counts in this table correspond to ours, the repeated instances for some measurement occasions (which could indicate data entry errors, but also be due to large variability in the time inteval between measurements) were not reported in the original analysis (i.e., Table 1). This suggests that the 2 participants with repeated occasions were simply counted and measured twice on one occasion and missing from another one.
Another good question to ask is:
- How does the number of
elapsed.dayscorrespond to the measurement occasions?
This brings us to a more general principle and corresponding practice that involves different ways of viewing distributions of values.
Viewing distributions
- Principle 6: Inspect the distributions of (continuous) variables.
Plotting the distribution of individual variables (e.g., with histograms or frequency polygons) provides a good overview over the values present in the data and any unusual or unexpected values.
Histograms
A histogram provides a cumulative overview of a variable’s values along one dimension (see our examples in visualizing data). Here, we use a histogram of elapsed.days to address the question:
- How are the measurement times (
elapsed.days) distributed overall (and relative to the stated times of eachoccasion)?
## Data:
df <- AHI_CESD
# From Table 1 (p. 4):
occ_days <- c(0, 7, 14, 38, 98, 189)
names(occ_days) <- c("0: pre-test", "1: post-test", "2: 1-week", "3: 2-weeks", "4: 3-months", "5: 6-months")
occ_days
#> 0: pre-test 1: post-test 2: 1-week 3: 2-weeks 4: 3-months
#> 0 7 14 38 98
#> 5: 6-months
#> 189
# (a) Histogram:
ggplot(df, aes(x = elapsed.days)) +
geom_histogram(fill = seeblau, binwidth = 1) +
geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) +
labs(title = "Distribution of occasions (a: histogram)") +
theme_bw()
Note: The first 3 occasions are as expected. However, occasions 4 to 6 appear shifted to the left (i.e., were about 7 days earlier than stated).
Alternative ways of viewing distributions
Alternative ways to plot cumulative distributions of values include frequency polygons, density plots, and rug plots.
# (b) frequency polygon:
ggplot(df, aes(x = elapsed.days)) +
# geom_histogram(fill = seeblau, binwidth = 1) +
geom_freqpoly(binwidth = 7, color = seeblau) +
geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) +
labs(title = "Distribution of occasions (b: frequency polygon)") +
theme_bw()
# (c) density plot:
ggplot(df, aes(x = elapsed.days)) +
# geom_histogram(fill = seeblau, binwidth = 1) +
geom_density(fill = seeblau) +
geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) +
labs(title = "Distribution of occasions (c: density plot)") +
theme_bw()
# (d) rug plot:
ggplot(df, aes(x = elapsed.days)) +
geom_freqpoly(binwidth = 7, color = seeblau) +
geom_rug(size = 1, color = "black", alpha = 1/4) +
geom_vline(xintercept = occ_days, color = "firebrick", linetype = 2) +
labs(title = "Distribution of occasions (d: frequency polygon with rug plot)") +
theme_bw()
Filtering values
Once we have detected something noteworthy or strange, we may want to mark or exclude some observations (rows) or variables (columns). As we can easily select rows and filter cases (see dplyr::select and dplyr::filter in the last session), it may be good to create and include some dedicated filter variables in our data. Let’s use the participant data posPsy_p_info to illustrate how we can create filter variables and then filter cases:
# Data:
p_info <- posPsy_p_info
dim(p_info) # 295 x 6
#> [1] 295 6In previous exercises (see Exercise 5 of WPA02 and Exercise 4 of WPA03), we answered the question:
- What is the
agerange of participants (overall and byintervention)?
by using both ggplot and dplyr pipes. For instance, we can answer questions about the distribution of age values by plotting histograms:
# Age range:
range(p_info$age)
#> [1] 18 83
summary(p_info$age)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 18.00 34.00 44.00 43.76 52.00 83.00
# Histogramm showing the overall distribution of age
ggplot(p_info) +
geom_histogram(mapping = aes(age), binwidth = 2, fill = "gold", col = "black") +
theme_bw() +
labs(title = "Distribution of age values (overall)")
# Create 4 histogramms showing the distribution of age by intervention:
ggplot(p_info) +
geom_histogram(mapping = aes(age), binwidth = 5, fill = seeblau, col = "black") +
theme_bw() +
labs(title = "Distribution of age values by intervention (all data)") +
facet_grid(.~intervention)
For practice purposes, suppose we only wanted to include participants up to an age of 70 years in some analysis. Rather than dropping these participants from the file, we can introduce a filter variable that is TRUE when some criterion (here: age > 70) is satisfied, and otherwise FALSE:
- Principle 7: Use filter variables to identify sub-sets of observations (rows).
# How many participants are over 70?
sum(p_info$age > 70) # 6 people with age > 70
#> [1] 6
# Which ones?
which(p_info$age > 70) # 51 83 114 155 215 244
#> [1] 51 83 114 155 215 244
# Show their details (using dplyr):
p_info %>%
filter(age > 70)
#> # A tibble: 6 x 6
#> id intervention sex age educ income
#> <int> <int> <int> <int> <int> <int>
#> 1 51 2 1 83 2 2
#> 2 83 4 1 71 2 1
#> 3 114 1 1 71 4 1
#> 4 155 3 2 71 4 1
#> 5 215 4 1 75 3 2
#> 6 244 2 1 77 3 2
# Add a corresponding filter variable to df:
p_info <- p_info %>%
mutate(over_70 = (age > 70))
dim(p_info) # => 7 variables (now including over_70 as a logical variable):
#> [1] 295 7
head(p_info)
#> # A tibble: 6 x 7
#> id intervention sex age educ income over_70
#> <int> <int> <int> <int> <int> <int> <lgl>
#> 1 1 4 2 35 5 3 FALSE
#> 2 2 1 1 59 1 1 FALSE
#> 3 3 4 1 51 4 3 FALSE
#> 4 4 3 1 50 5 2 FALSE
#> 5 5 2 2 58 5 2 FALSE
#> 6 6 1 1 31 5 1 FALSE
# Check details again (but applying filter to over_70):
p_info %>%
filter(over_70)
#> # A tibble: 6 x 7
#> id intervention sex age educ income over_70
#> <int> <int> <int> <int> <int> <int> <lgl>
#> 1 51 2 1 83 2 2 TRUE
#> 2 83 4 1 71 2 1 TRUE
#> 3 114 1 1 71 4 1 TRUE
#> 4 155 3 2 71 4 1 TRUE
#> 5 215 4 1 75 3 2 TRUE
#> 6 244 2 1 77 3 2 TRUEIn the present case, filter(over_70) is about as long and complicated as filter(age > 70) and thus not really necessary. However, defining explicit filter variables can pay off when constructing more complex filter conditions or when needing several sub-sets of the data (e.g., for cross-validation purposes). Given an explicit filter variable, we can later filter any analysis (or plot) on the fly:
# Age distribution of participants up to 70 (i.e., not over 70):
p_info %>%
filter(over_70 == FALSE) %>%
ggplot() +
geom_histogram(mapping = aes(age), binwidth = 5, fill = "olivedrab3", col = "black") +
theme_bw() +
labs(title = "Distribution of age values by intervention (without participants over 70)") +
facet_grid(.~intervention)
# Alternatively, we can quickly create sub-sets of the data:
p_info_young <- p_info %>%
filter(over_70 == FALSE)
p_info_old <- p_info %>%
filter(over_70 == TRUE)
dim(p_info_young) # 289 participants (295 - 6) x 7 variables
#> [1] 289 7
dim(p_info_old) # 6 participants x 7 variables
#> [1] 6 7Viewing relationships
Most research hypotheses involve relationships or measures of correspondence between 2 or more variables. This brings us to another principle:
- Principle 8: Inspect relationships between (continuous or categorical) variables.
Scatterplots
A scatterplot shows the relationship between 2 (typically continuous) variables:
# Data:
df <- AHI_CESD
dim(df) # 992 50
#> [1] 992 50
# df
# Scatterplot (overall):
ggplot(df) +
geom_point(aes(x = ahiTotal, y = cesdTotal), size = 2, alpha = 1/4) +
geom_abline(intercept = 100, slope = -1, col = seeblau) +
labs(title = "Relationship between ahiTotal and cesdTotal (overall)",
x = "ahiTotal", y = "cesdTotal") +
theme_bw()
# Scatterplot (with facets):
ggplot(df) +
geom_point(aes(x = ahiTotal, y = cesdTotal), size = 1, alpha = 1/3) +
geom_abline(intercept = 100, slope = -1, col = seeblau) +
labs(title = "Relationship between ahiTotal and cesdTotal (by intervention)",
x = "ahiTotal", y = "cesdTotal") +
facet_wrap(~intervention) +
theme_bw()
Bar plots and line graphs
A relationship often of interest in psychology and other sciences are trends or developments over time. When studying trends or developments, the time variable is not always explicitly expressed in units of time (e.g., days, weeks, months), but often encoded implicitly (e.g., as repeated measurements).
- Principle 9: Inspect trends over time or repeated measurements.
Trends can be expressed in many ways. Any data that measures some variable more than once and uses some unit of time (e.g., in AHI_CESD we have a choice between occasion vs. elapsed.days) can be used to address the question:
- How do values change over time?
To visualize trends, we typically map the temporal variable (in units of time or the counts of repeated measurements) to the x-axis of a plot. If this is the case, many different types of plots can express developments over time. In our session on in visualizing data, we have already encountered bar plots and line graphs, which are well-suited to plot trends. Actually, our bar plots above (using either the raw data of AHI_CESD or our summary table id_by_occ to show the number of participants by occasion) expressed a trend: How many participants dropped out of the study from occasion 0 to 6. We can make the same point in a line plot:
# Data:
df <- AHI_CESD
dim(df) # 992 50
#> [1] 992 50
df
#> # A tibble: 992 x 50
#> id occasion elapsed.days intervention ahi01 ahi02 ahi03 ahi04 ahi05
#> <int> <int> <dbl> <int> <int> <int> <int> <int> <int>
#> 1 1 0 0 4 2 3 2 3 3
#> 2 1 1 11.8 4 3 3 4 3 3
#> 3 2 0 0 1 3 4 3 4 2
#> 4 2 1 8.02 1 3 4 4 4 3
#> 5 2 2 14.3 1 3 4 4 4 3
#> 6 2 3 32.0 1 3 4 4 4 4
#> 7 2 4 92.2 1 3 3 2 3 3
#> 8 2 5 182. 1 3 3 3 4 2
#> 9 3 0 0 4 3 3 2 4 2
#> 10 3 2 16.4 4 3 3 3 4 4
#> # ... with 982 more rows, and 41 more variables: ahi06 <int>, ahi07 <int>,
#> # ahi08 <int>, ahi09 <int>, ahi10 <int>, ahi11 <int>, ahi12 <int>,
#> # ahi13 <int>, ahi14 <int>, ahi15 <int>, ahi16 <int>, ahi17 <int>,
#> # ahi18 <int>, ahi19 <int>, ahi20 <int>, ahi21 <int>, ahi22 <int>,
#> # ahi23 <int>, ahi24 <int>, cesd01 <int>, cesd02 <int>, cesd03 <int>,
#> # cesd04 <int>, cesd05 <int>, cesd06 <int>, cesd07 <int>, cesd08 <int>,
#> # cesd09 <int>, cesd10 <int>, cesd11 <int>, cesd12 <int>, cesd13 <int>,
#> # cesd14 <int>, cesd15 <int>, cesd16 <int>, cesd17 <int>, cesd18 <int>,
#> # cesd19 <int>, cesd20 <int>, ahiTotal <int>, cesdTotal <int>
knitr::kable(id_by_occ)| occasion | n |
|---|---|
| 0 | 295 |
| 1 | 147 |
| 2 | 157 |
| 3 | 139 |
| 4 | 134 |
| 5 | 120 |
# Plot count of IDs by occasion as line graph:
# Line plot: Using id_by_occ (from above) as input:
line_dropouts <- ggplot(id_by_occ, aes(x = occasion)) +
geom_line(aes(y = n), size = 1.5, color = "gold") +
geom_point(aes(y = n), shape = 21, size = 4, stroke = 1, color = "black", fill = "gold") +
labs(title = "Number of participants by occasion (line)") +
theme_bw()
line_dropouts
When comparing this plot with our corresponding bar plot above, we note that ggplot has automatically scaled the y-axis to a range around our frequency values n (here: from about 110 to 300). As a consequence, the number of dropouts over occasions looks more dramatic in the line graph than it did in the bar plot. We can correct for this difference by explicitly scaling the y-axis of the line graph or by plotting both geoms in one plot. And whenever combining multiple geoms in the same plot, we can move common aesthetic elements (e.g., y = n) to the 1st line of our ggplot command, so that we don’t have to repeat them and can only list the specific aesthetics of each geom later:
## Data:
# knitr::kable(id_by_occ)
# Line plot (from above) with corrected y-scale:
line_dropouts +
scale_y_continuous(limits = c(0, 300))
# Multiple geoms:
# Bar + line plot: Using id_by_occ (from above) as input:
ggplot(id_by_occ, aes(x = occasion, y = n)) +
geom_bar(stat = "identity", fill = seeblau) +
geom_line(size = 1.5, color = "gold") +
geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "gold") +
labs(title = "Number of participants by occasion (bar + line)") +
theme_bw()
Two key questions in the context of our study are:
- How do participants’ scores of happiness (
ahiTotal) and depression (cesdTotal) change over time?
- Do these trends over time vary by a participant’s type of
intervention?
Visualizing these trends will take us more than half-way towards our conclusion, except for a few statistical tests. To answer the first question for the happiness scores (ahiTotal), we compute our desired values with dplyr and then plot the resulting table as a line graph (with an automatic and a full y-scale):
# Data:
df <- AHI_CESD
dim(df) # 992 50
#> [1] 992 50
# df
# Table 1:
# Frequency n, mean ahiTotal, and mean cesdTotal by occasion:
tab_scores_by_occ <- df %>%
group_by(occasion) %>%
summarise(n = n(),
ahiTotal_mn = mean(ahiTotal),
cesdTotal_mn = mean(cesdTotal)
)
knitr::kable(head(tab_scores_by_occ))| occasion | n | ahiTotal_mn | cesdTotal_mn |
|---|---|---|---|
| 0 | 295 | 69.70508 | 15.06441 |
| 1 | 147 | 72.35374 | 12.96599 |
| 2 | 157 | 73.04459 | 12.36943 |
| 3 | 139 | 74.55396 | 11.74101 |
| 4 | 134 | 75.19403 | 11.70896 |
| 5 | 120 | 75.83333 | 12.83333 |
# Line graph of scores over occasions:
plot_ahi_trend <- ggplot(tab_scores_by_occ, aes(x = occasion, y = ahiTotal_mn)) +
# geom_bar(aes(y = n), stat = "identity", fill = seeblau) +
geom_line(size = 1.5, color = "forestgreen") +
geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "forestgreen") +
labs(title = "ahiTotal_mn by occasion (line)") +
theme_bw()
plot_ahi_trend
# Line graph with corrected y-scale:
plot_ahi_trend +
scale_y_continuous(limits = c(0, 80))
Practice: Create the same plots for mean depression scores cesdTotal. Can you plot both lines in the same plot? What do you see?
To study the trends of scores by participant’s type of intervention, we need to slightly adjust our workflow by adding the intervention variable into our summary table, so that we can later use it in our plotting commands. As before, we will take care of happiness (ahiTotal) and leave it to you to practice on depression (cesdTotal):
# Data:
df <- AHI_CESD
dim(df) # 992 50
#> [1] 992 50
# df
# Table 2:
# Scores (n, ahiTotal_mn, and cesdTotal_mn) by occasion AND intervention:
tab_scores_by_occ_iv <- df %>%
group_by(occasion, intervention) %>%
summarise(n = n(),
ahiTotal_mn = mean(ahiTotal),
cesdTotal_mn = mean(cesdTotal)
)
dim(tab_scores_by_occ_iv) # 24 x 5
#> [1] 24 5
knitr::kable(head(tab_scores_by_occ_iv)) # print table| occasion | intervention | n | ahiTotal_mn | cesdTotal_mn |
|---|---|---|---|---|
| 0 | 1 | 72 | 68.38889 | 15.05556 |
| 0 | 2 | 76 | 68.75000 | 16.21053 |
| 0 | 3 | 74 | 70.18919 | 16.08108 |
| 0 | 4 | 73 | 71.50685 | 12.84932 |
| 1 | 1 | 30 | 69.53333 | 15.30000 |
| 1 | 2 | 48 | 71.58333 | 14.58333 |
# Line graph of scores over occasions AND interventions:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = ahiTotal_mn,
group = intervention)) +
geom_line(size = 1.5, color = "forestgreen") +
geom_point(shape = 21, size = 4, stroke = 1, color = "black", fill = "forestgreen") +
labs(title = "ahiTotal_mn by occasion and intervention (a: lines + points)") +
theme_bw()
Note that we added group = intervention to the general aesthetics to instruct geom_line to interpret instances of the same intervention but different occasions as belonging together (i.e., to the same line). When using multiple lines, our previous color settings (in geom_line and geom_point) no longer make sense. Instead, we would like to vary the color of our lines and points by intervention. We can achieve this by also moving these aesthetics to the general aes() of the plot:
# Line graph of scores over occasions AND interventions:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = ahiTotal_mn,
group = intervention,
color = intervention,
fill = intervention)) +
geom_line(size = 1.5) +
geom_point(shape = 21, size = 4, stroke = 1, color = "black") +
labs(title = "ahiTotal_mn by occasion and intervention (b: lines + points)") +
theme_bw()
The resulting plot is still suboptimal. The automatic choice of a color scale (in “shades of blue”) and the corresponding legend (on the right) indicates that intervention (i.e., a variable of type integer in df and tab_scores_by_occ_iv) was interpreted as a continuous variable. To turn it into a discrete variable, we can use factor(intervention) for both the color and the fill aesthetics:
# Line graph of scores over occasions AND interventions:
plot_ahi_by_occ_iv <- ggplot(tab_scores_by_occ_iv,
aes(x = occasion, y = ahiTotal_mn,
group = intervention,
color = factor(intervention),
fill = factor(intervention))) +
geom_line(size = 1.5) +
geom_point(shape = 21, size = 4, stroke = 1, color = "black") +
labs(title = "ahiTotal_mn by occasion and intervention (c: lines + points, IVs by color)") +
theme_bw()
plot_ahi_by_occ_iv
An more honest version of the same plot could manually set a different color scheme, but would still look a lot more sobering:
# Line graph of scores over occasions AND interventions:
plot_ahi_by_occ_iv +
scale_y_continuous(limits = c(0, 80)) +
scale_color_brewer(palette = "Set1") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Mean ahiTotal scores by occasion and intervention (d)",
x = "Occasion", y = "Happiness (mean ahiTotal)",
color = "Intervention:", fill = "Intervention:")
Practice: Create the same plots for mean depression scores cesdTotal. What do you see?
When viewing and contrasting the trends of groups, we can always zoom in further:
- How does this look for individual participants?
Well, how nice of you to ask – so let’s see:
# Data:
df <- AHI_CESD
dim(df) # 992 50
#> [1] 992 50
# df
# +++ here now +++Jitter, box, and violin plots
A measure of correspondence that is common in psychology asks whether the values of some continuous variable vary as a function of the levels of a categorical variable. With regard to our data in AHI_CESD we may wonder:
- Do the happiness scores (
ahiTotal) vary byintervention?
To visualize the relationship, we cannot use a scatterplot with the mapping x = intervention and y = ahiTotal, as there would only be 4 distinct values for x (go ahead plotting it, if you want to see it). Fortunately, there’s a range of alternatives that allow plotting the raw values and distributions of a continuous variable as a function of a categorical one:
## Data:
# df <- AHI_CESD
# dim(df) # 992 50
# df
# (a) Jitterplot:
ggplot(df) +
geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) +
labs(title = "Values of ahiTotal by intervention (a: jitter)",
x = "intervention", y = "ahiTotal") +
theme_bw()
# (b) Box plot:
ggplot(df) +
geom_boxplot(aes(x = factor(intervention), y = ahiTotal), color = seeblau, fill = "grey95") +
labs(title = "Values of ahiTotal by intervention (b: boxplot)",
x = "intervention", y = "ahiTotal") +
theme_bw()
# Note that we use factor(intervention), as intervention is an integer (i.e., contiuous) variable.
# (c) Violin plot:
ggplot(df) +
geom_violin(aes(x = factor(intervention), y = ahiTotal), size = 1.5, color = seeblau, fill = "whitesmoke") +
labs(title = "Values of ahiTotal by intervention (c: violin)",
x = "intervention", y = "ahiTotal") +
theme_bw()
Note that we sometimes use factor(intervention), as intervention in is an integer (i.e., contiuous) variable.
Sometimes combining 2 or more geoms can be more informative than just using one:
# (d) Combining jitter with boxplot:
ggplot(df) +
geom_boxplot(aes(x = factor(intervention), y = ahiTotal), color = seeblau, fill = "grey95", alpha = 1) +
geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) +
labs(title = "Values of ahiTotal by intervention (d: jittered boxes)",
x = "intervention", y = "ahiTotal") +
theme_bw()
# (e) Combining violin plot with jitter:
ggplot(df) +
geom_violin(aes(x = factor(intervention), y = ahiTotal), size = 1.5, color = seeblau) +
geom_jitter(aes(x = intervention, y = ahiTotal), width = .1, size = 2, alpha = 1/4) +
labs(title = "Values of ahiTotal by intervention (e: jittered violins)",
x = "intervention", y = "ahiTotal") +
theme_bw()
# (f) Combining violin plot with boxplot and jitter:
ggplot(df, aes(x = factor(intervention), y = ahiTotal, color = factor(intervention))) +
geom_violin(size = 1.5, color = seeblau) +
geom_boxplot(width = .30, color = "grey20", fill = "grey90") +
geom_jitter(width = .05, size = 2, alpha = 1/3) +
labs(title = "Values of ahiTotal by intervention (f: jittered violins with boxes)",
x = "intervention", y = "ahiTotal", color = "Intervention:") +
scale_color_brewer(palette = "Set1") +
theme_bw()
When using multiple geoms in one plot, their order matters, as later geoms are printed on top of earlier ones. In addition, combining geoms typically requires playing with aesthetics. Incidentally, note how the last plot moved some redundant aesthetic mappings (i.e., x, y, and color) from the geoms to the 1st line of the command (i.e., from the mapping argument of the geoms to the general mapping argument).
Tile plots
We already discovered that bar plots can either count cases or show some pre-computed values (when using stat = "identity", see Bar plots in Visualizing data. In the following, we show that we can also compute summary tables and then display their values by the color fill gradient of tiles, or by the size of points:
# Frequency of observations by occasion (x) and intervention (y):
# Data:
df <- AHI_CESD
# dim(df) # 992 50
# df
# Table 1:
# Frequency n, mean ahiTotal, and mean cesdTotal by occasion:
# Use tab_scores_by_occ (from above):
knitr::kable(head(tab_scores_by_occ))| occasion | n | ahiTotal_mn | cesdTotal_mn |
|---|---|---|---|
| 0 | 295 | 69.70508 | 15.06441 |
| 1 | 147 | 72.35374 | 12.96599 |
| 2 | 157 | 73.04459 | 12.36943 |
| 3 | 139 | 74.55396 | 11.74101 |
| 4 | 134 | 75.19403 | 11.70896 |
| 5 | 120 | 75.83333 | 12.83333 |
# Table 2:
# Scores (n, ahiTotal_mn, and cesdTotal_mn) by occasion AND intervention:
# Use tab_scores_by_occ_iv (from above):
knitr::kable(head(tab_scores_by_occ_iv))| occasion | intervention | n | ahiTotal_mn | cesdTotal_mn |
|---|---|---|---|---|
| 0 | 1 | 72 | 68.38889 | 15.05556 |
| 0 | 2 | 76 | 68.75000 | 16.21053 |
| 0 | 3 | 74 | 70.18919 | 16.08108 |
| 0 | 4 | 73 | 71.50685 | 12.84932 |
| 1 | 1 | 30 | 69.53333 | 15.30000 |
| 1 | 2 | 48 | 71.58333 | 14.58333 |
# Tile plot 1:
# Frequency (n) of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
geom_tile(aes(fill = n), color = "white", size = .1) +
geom_text(aes(label = n), color = "white") +
labs(title = "Frequency of each intervention by occasion (tile)",
x = "Occasion", y = "Intervention", fill = "Frequency:") +
theme_bw()
# Tile plot 2:
# ahiTotal_mn of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
geom_tile(aes(fill = ahiTotal_mn), color = "grey75", size = 1) +
labs(title = "Mean ahiTotal score of each intervention by occasion (tile)",
x = "Occasion", y = "Intervention", fill = "ahiTotal_mn:") +
scale_fill_gradient(low = "black", high = "gold") +
theme_bw()
Note the fancy color gradient of the last tile plot: The lowest average scores of happiness (or ahiTotal_mn) are shown in "black", the highest ones in "gold". As the scores vary on a continuous scale, using a continuous color scale is appropriate here. And assigning high happiness to “gold” hopefully makes the plot easier to interpret.
Practice: Create a tile plot that shows the developments of mean depression scores cesdTotal by measurement occasion and intervention. Choose a color scale that you find appropriate. What do you see?
Hint: Your result could look like the following plot:
#> # A tibble: 24 x 5
#> # Groups: occasion [6]
#> occasion intervention n ahiTotal_mn cesdTotal_mn
#> <int> <int> <int> <dbl> <dbl>
#> 1 0 1 72 68.4 15.1
#> 2 0 2 76 68.8 16.2
#> 3 0 3 74 70.2 16.1
#> 4 0 4 73 71.5 12.8
#> 5 1 1 30 69.5 15.3
#> 6 1 2 48 71.6 14.6
#> 7 1 3 24 74 12.6
#> 8 1 4 45 74.2 9.87
#> 9 2 1 38 70.3 13.6
#> 10 2 2 48 73.0 12.2
#> # ... with 14 more rows
Point size plots
Another way of expressing the value of a continuous variable – like a frequency n or an average score ahiTotal_mn – as a function of 2 categorical variables is by mapping it to the size dimension of a point. In the following, we illustrate this by creating the same plot twice: First from raw data (df, using geom_count) and then from our summary table (tab_scores_by_occ_iv from above, using geom_point with size mapped to our frequency count n):
# Point plot 1a:
# Count frequency of values for each interventin by occasion:
ggplot(df, aes(x = occasion, y = intervention)) +
geom_count(color = seeblau) +
labs(title = "Frequency of each intervention by occasion (count)",
x = "Occasion", y = "Intervention", size = "Count:") +
theme_bw()
# Point plot 1b:
# Frequency (n mapped to size) of each interventin by occasion:
ggplot(tab_scores_by_occ_iv, aes(x = occasion, y = intervention)) +
geom_point(aes(size = n), color = seeblau) +
labs(title = "Frequency of each intervention by occasion (point size)",
x = "Occasion", y = "Intervention", size = "Frequency:") +
theme_bw()
Mixing geoms with multiple aesthetics
By combining the geom_tile and geom_point plots from above with geom_text(aes(label = n), ... we could simultaneously express 2 different continuous variables on 2 different dimensions (here fill color or size vs. the value shown by a label of text):
# Combo plot 1:
# ahiTotal_mn of each interventin by occasion:
ggplot(summary, aes(x = occasion, y = intervention)) +
geom_tile(aes(fill = ahiTotal_mn), color = "grey95", size = 1) +
geom_text(aes(label = n), color = "white", size = 5, fontface = 2) +
labs(title = "Mean ahiTotal score of each intervention by occasion (tile)",
x = "Occasion", y = "Intervention", fill = "ahiTotal_mn:") +
scale_fill_gradient(low = "grey30", high = "gold") +
theme_bw()
# Combo plot 2b:
# ahiTotal_mn (mapped to point size) and n (text label) for each intervention x occasion:
ggplot(summary, aes(x = occasion, y = intervention)) +
geom_point(aes(size = ahiTotal_mn), color = "forestgreen") +
geom_text(aes(label = n), color = "white", size = 3, fontface = 2) +
labs(title = "ahiTotal_mn and n by intervention x occasion (point size + text label)",
x = "Occasion", y = "Intervention", size = "Frequency:") +
scale_size_continuous(range = c(5, 11)) +
theme_bw()
However, this mix of multiple dimensions tends to get too confusing, even for ggplot enthusiasts. So let’s conclude by a cautionary note.
Plotting with a purpose
Our last examples have shown that we should better create 2 separate graphics when making multiple points or expressing more than a single relationship in a plot. This brings us to:
- Principle 10: A plot should convey its message as clearly as possible.
So go ahead and use the awesome ggplot machine to combine multiple geoms and adjust their aesthetics for as long as this makes your plot prettier. However, rather than getting carried away by the plethora of options, always keep in mind the message that you want to convey. If additional fiddling with aesthetics does not help to clarify your point, you are probably wasting your time by trying to do too much.
Exercises (WPA04)
Working through the above examples should have provided you with a pretty clear picture of the data in AHI_CESD. However, we haven’t touched the data contained in posPsy_wide yet. Given that this is described as a modified and corrected version of AHI_CESD, we should not expect to find completely different results in it. So let’s use posPsy_wide to exercise and improve our skills in EDA.
More on EDA
- study
vignette("ggplot")and the documentation forggplotand various geoms (e.g.,geom_); - study https://ggplot2.tidyverse.org/reference/ and its examples;
- see the cheat sheet on data visualization;
- read Chapter 3: Data visualization and Chapter 7: Exploratory data analysis (EDA) and complete their exercises.
Conclusion
All ds4psy essentials:
| Nr. | Topic |
|---|---|
| 0. | Syllabus |
| 1. | Basic R concepts and commands |
| 2. | Visualizing data |
| 3. | Transforming data |
| 4. | Exploring data (EDA) |
| +. | Datasets |
[Last update on 2018-11-28 11:48:12 by hn.]
The participant part of this data has been used as
p_infoin the previous sessions.↩Actually, listing everything required to execute a file can get quite excessive (check out
sessionInfo()or thesession_infoandpackage_infofunctions ofdevtools, in case you’re interested, and consider using thepackratpackage in case you want to preserve or share your current setup). Hence, most people only explicitly list and load non-standard packages (in their current version).↩